      subroutine blusol(ig)

      use param
      use allocs

#include <global.h>
      include 'geom.h'
      include 'setls.h'
      include 'eigv.h'

      integer,intent(in) :: ig
      integer tn
      real(NBF) invd
      real(NBF),dimension(2) :: bmat

!     forward substitution      
!     define pseudo zeroth flux
      psphi(1:2,1,ig)=src(1:2,1,ig) 

      do j=2,tnode
        ja=j-1
        bmat(1:2)=src(1:2,j,ig)
        psphi(1,j,ig)=bmat(1)-bdl(1,j,ig)*psphi(1,ja,ig)- &   
                     bdl(2,j,ig)*psphi(2,ja,ig)
       psphi(2,j,ig)=bmat(2)-bdl(3,j,ig)*psphi(1,ja,ig)- &
                     bdl(4,j,ig)*psphi(2,ja,ig)
      enddo

!     backward substitution
      bmat(1:2)=psphi(1:2,tnode,ig)
      invd=1/(bdd(1,tnode,ig)*bdd(4,tnode,ig)-bdd(2,tnode,ig)*bdd(3,tnode,ig))
      psphi(1,tnode,ig)=invd*(bdd(4,tnode,ig)*bmat(1)-bdd(2,tnode,ig)*bmat(2))
      psphi(2,tnode,ig)=invd*(bdd(1,tnode,ig)*bmat(2)-bdd(3,tnode,ig)*bmat(1))
      do j=tnode-1,1,-1
        jb=j+1
        bmat(1:2)=psphi(1:2,j,ig)
        bmat(1)=bmat(1)-bdu(1,j,ig)*psphi(1,jb,ig)- &
               bdu(2,j,ig)*psphi(2,jb,ig)
        bmat(2)=bmat(2)-bdu(3,j,ig)*psphi(1,jb,ig)- &
               bdu(4,j,ig)*psphi(2,jb,ig)
        invd=1/(bdd(1,j,ig)*bdd(4,j,ig)-bdd(2,j,ig)*bdd(3,j,ig))
        psphi(1,j,ig)=invd*(bdd(4,j,ig)*bmat(1)-bdd(2,j,ig)*bmat(2))
        psphi(2,j,ig)=invd*(bdd(1,j,ig)*bmat(2)-bdd(3,j,ig)*bmat(1))
      enddo

      end subroutine
